% Figure O6
axoptions={'scaled ticks = false',...
           'y tick label style={/pgf/number format/.cd, fixed, fixed zerofill,precision=0, set thousands separator={}}',...
           'x tick label style={/pgf/number format/.cd,precision=1, set thousands separator={}}',...
           'legend style={font=\normalsize}'}; %\scriptsize
load ../data/DataforMatlab/IP.csv
[UB,LB] = Cal_B(I_vec, B_vec, pvec);
%% Figure O6:
wT = I_vec(:,:,T);
wT(isnan(wT)) = [];
[wT, index] = unique(wT);
%%Percentile
Wfull = reshape(I_vec(:,:,:),[N,T]);
U = reshape(U_vec,[N,T]);

tmp = sum([Wfull <=IP(1,:)]);
tmp(tmp ==0 ) = 1;
PC(1,:) = tmp;
PC(2,:) = sum([Wfull <=IP(2,:)]);
PC(3,:) = sum([Wfull <=IP(3,:)]);
PC(4,:) = sum([Wfull <=IP(4,:)]);
PC(5,:) = sum([Wfull <=IP(5,:)]);
PC(6,:) = sum([Wfull <=IP(6,:)]);
PC(7,:) = sum([Wfull <=IP(7,:)]);
PC(8,:) = sum([Wfull <=IP(8,:)]);
PC(9,:) = sum([Wfull <=IP(9,:)]);
PC(PC==0)=1;
Wq = zeros(9,T);
Uq = zeros(9,T);
for t = 1:T
    for p = 1:9
        Wq(p,t) = Wfull(PC(p,t),t);
        Uq(p,t) = U(PC(p,t),t);
    end
end

U = reshape(U_vec,[N,T]);

[~, indmax] = max(U);
[~, indmin] = min(U);
W = Wfull;

for tau = 1:T
Wmax(1,tau) = Wfull(indmax(tau),tau);
Wmin(1,tau) = Wfull(indmin(tau),tau);
wind = (indmin(tau):indmax(tau));
W(:,tau) = NaN;
W(wind,tau) = Wfull(wind,tau);
end
logUmax = max(log(U));
logUmin = min(log(U));

OurxaxisUK = [Wmin(:,end);Wq(:,end);Wmax(:,end)]*52;
OuryaxisUK = [exp(logUmin(:,end)+log(52)) ;(Uq(:,end)*52);exp(logUmax(:,end)+log(52)) ];


UBscale = griddedInterpolant(wT*52,UB(index)*52,'linear','none');
LBscale = griddedInterpolant(wT*52,LB(index)*52,'linear','none');
Aggfun = griddedInterpolant(wT*52,LB(index)*52,'linear','none');
figUB = UBscale(OurxaxisUK);
figLB = LBscale(OurxaxisUK);
figure('name','Figure A.6','NumberTitle','off')
loglog(OurxaxisUK,OuryaxisUK,'LineWidth', 2); hold on
loglog(OurxaxisUK,figUB,'--','LineWidth', 1); 
loglog(OurxaxisUK,figLB,'--','LineWidth', 1); 
ax.XAxis.Exponent = 0;
ax.XAxis.TickLabelFormat = '%.0f';
xlim([min(OurxaxisUK), 100000])
ylim([min(figLB), max(figUB)])
xticklabels({'10,000','100,000'})
yticklabels({'1,000','10,000'})


xlabel('Income in 2017')
ylabel('1974 Income that gives the same utility as in 2017')
legend('Money metric','Upper Bound','Lower Bound','Location','northwest')
%matlab2tikz('../fig/UKBrandell.tex','extraAxisOptions',axoptions)



function [UB,LB] = Cal_B(I_vec, b_vec, pvec)

T = size(b_vec,3);N = size(b_vec,2);I = size(b_vec,1);
% To obtain Money-Metrics instead of cost-of-living, the algorithm is
% applied from T to 1. 
b_vec = flip(b_vec,3);
pvec = flip(pvec,3);
I_vec = flip(I_vec,3);

pq = I_vec.*b_vec;
q  = pq./pvec;

q1 = q(:,:,1);
ptq1 = sum(q1.*pvec,1);
bound = 1/100;
diff = 10;
qis = zeros(I,N,T);
qis(:,:,1) = q(:,:,1);
Fs = q;
count = 0;
% Algorithm A

for t = 2:T
    wt = I_vec(:,:,t);
    wt(isnan(wt)) = [];
    qt = q(:,:,t);
    qt = qt(:,1:size(wt,2));
    [wt, index] = unique(wt);

    for i = 1:I
        qit_tmp = qt(i,:);
        qit = griddedInterpolant(wt,qit_tmp(index),'linear','none');
        qis(i,:,t) = qit(ptq1(:,:,t));
    end
end
Fs = qis;
while diff > bound 
% updating

for t = 2:T
    for i = 1:I
        qit_tmp = qt(i,:);
        qit = griddedInterpolant(wt,qit_tmp(index),'linear','none');
        qis(i,:,t) = qit(min(sum(Fs.*pvec(:,:,t),1),[],3));
    end
end
Fs_1  = qis;
err = abs(Fs - Fs_1);
diff = max(max(max(err(:,1:N,:),[],3)));
count = count + 1;
Fs = Fs_1;

end
% upper bound
UB = min(sum(Fs_1(:,:,2:end).*pvec(:,:,T),1),[],3);


%% B
q  = pq./pvec;
q1 = q(:,:,1);
ptq1 = sum(q1.*pvec,1);
bound = 1/100;
diff = 10;
qis = zeros(I,N,T);
qis(:,:,1) = q(:,:,1);
Fs = q;
%v1 = v_vec(:,:,1);
count = 0;
dummy0 = zeros(1,N);
% step 1
for t = 2:T
    pq_b = sum(pvec(:,:,1).*q(:,:,t),1);
    pq_b(isnan(pq_b)) = [];
    wt = I_vec(:,:,t);
    wt = wt(:,1:size(pq_b,2));
    qt = q(:,:,t);
    qt = qt(:,1:size(pq_b,2));
    [pq_b, index] = unique(pq_b);
    p1qt_inv = griddedInterpolant(pq_b,wt(index),'linear','none');
    x = p1qt_inv(sum(pvec(:,:,1).*q(:,:,1),1));

    for i = 1:I
        qit_tmp = qt(i,:);
        qit_tmp2 = qit_tmp(index);
        wt2 = wt(index);
        [wt2, index2] = unique(wt2);
        qit = griddedInterpolant(wt2,qit_tmp2(index2),'linear','none');
        qis(i,:,t) = qit(x);
    end

end
% updating
Fs  = qis;
% step 2
while diff > bound
for t = 2:T
    xs = zeros(1,N,T);
for s = 2:T
        pqt_b = sum(pvec(:,:,s).*q(:,:,t),1);
        pqt_b(isnan(pqt_b)) = [];
        wt = I_vec(:,:,t);
        wt = wt(:,1:size(pqt_b,2));
        [pqt_b2, index2] = unique(pqt_b);
            %%%
        psqt_inv = griddedInterpolant(pqt_b2,wt(index2),'linear','none');
        xs(:,:,s) = psqt_inv(sum(pvec(:,:,s).*Fs(:,:,s),1));
end
x = max(xs(:,:,s),[],3);
    for i = 1:I
            wt = I_vec(:,:,t);
            wt(isnan(wt)) = [];
            qt = q(:,:,t);
            qt = qt(:,1:size(wt,2));
            [wt, index] = unique(wt);
        qit_tmp = qt(i,:);
        qit_tmp2 = qit_tmp(index);
        qit = griddedInterpolant(wt,qit_tmp2,'linear','none');
        qis(i,:,t) = qit(x);
    end
end
Fs_1  = qis;
err = abs(Fs - Fs_1);
diff = max(max(max(err(:,1:N,:),[],3)));
Fs = Fs_1;
count = count + 1;
end
LB = min(sum(Fs_1(:,:,2:end).*pvec(:,:,T),1),[],3);

end

